/* 
 * File:   inicial.h
 * Author: jonk
 *
 * Created on 15 de enero de 2012, 20:50
 */

#ifndef INICIAL_H
#define INICIAL_H

void SetRandomPart(particula *part, parametros *param)
{
    int nok=0;
    bool good;
    double r2;
    double dx,dy,dz;

    #ifdef POLYDISPERSE
        double poly=0.3;
        double Vpart=0.0;
        double u1,u2;
        ///Calculamos el tamanyo de caja para la polidispersidad considerada
        for(int i=0;i< param->npart;i++)
        {
            u1=double(rand())/RAND_MAX;
            u2=double(rand())/RAND_MAX;
            //Aqui introduciremos la distribucion de tamaños
            part[i].diameter=exp(poly*sqrt(-2.*log(u1))*cos(2.*pi*u2));
            Vpart+=pi*pow(part[i].diameter,3.0)/6.0;
            part[i].mass=pow(part[i].diameter,3.0);
        }
        param->vol=Vpart/phi;
        #ifdef CONFINED
            if (param->h0>0.)
            {
                param->lbox[0]=pow(param->vol/h0,1./3.);
                param->lbox[2]=param->h0*param->lbox[0];
                param->lbox[1]=param->lbox[0];
            }
            else{
                param->lbox[0]=pow(param->vol,1./3.);
                param->lbox[1]=param->lbox[2]=param->lbox[0];
            }
        #endif
        #ifndef CONFINED
            param->lbox[0]=pow(param->vol,1./3.);
            param->lbox[1]=param->lbox[2]=param->lbox[0];
        #endif
#else 
    double Vpart=0.0;
    for(int i=0;i< param->npart;i++)
    {
        part[i].diameter=1.0;
        Vpart+=pi*pow(part[i].diameter,3.0)/6.0;
        part[i].mass=pow(part[i].diameter,3.0);
    }
    param->vol=Vpart/param->phi;
    #ifdef CONFINED
        if (param->h0>0.)
        {
            param->lbox[0]=pow(param->vol/param->h0,1./3.);
            param->lbox[2]=param->h0*param->lbox[0];
            param->lbox[1]=param->lbox[0];
        }
        else
        {
            param->lbox[0]=pow(param->vol,1./3.);
            param->lbox[1]=param->lbox[2]=param->lbox[0];
//            cout<<param->lbox[0]<<endl;
        }
    #else
        param->lbox[0]=pow(param->vol,1./3.);
        param->lbox[1]=param->lbox[2]=param->lbox[0];
    #endif
#endif

    for(int i=0;i< param->npart;i++)
    {
        part[i].pos[0]=param->lbox[0]*(double(rand())/RAND_MAX);
        part[i].pos[1]=param->lbox[1]*(double(rand())/RAND_MAX);
        part[i].pos[2]=(param->lbox[2]-2.*part[i].diameter)*(double(rand())/RAND_MAX)+part[i].diameter;

        BC(&part[i],*param);
        ///Comprobamos solapamiento con otras particulas
        good=true;
        for(int j=0;j<nok;j++)
        {
            ///Recorremos todas las particulas que ya esten metidas en la lista
            if(good)
            {
                r2=distancia(part[i],part[j], &dx, &dy, &dz, *param);
                #ifdef POLYDISPERSE
                    if(sqrt(r2)<0.5*(part[i].radio+part[j].radio))
                    {
                        ///Hay solapamiento entre las particulas i y j. No aceptamos la particula i
                        good=false;
                        i--;
                    }
                #else
                    if(r2<1.2){
                        ///Hay solapamiento entre las particulas i y j. No aceptamos la particula i
                        good=false;
                        i--;
                    }
                #endif
                }
            }
            if(good)
            {
                ///Si pasa por aqui la particula no se solapa con ninguna otra --> se acepta
                cout<<"Numero de particulas colocadas: "<<nok<<endl;
                nok++;
            }
    }
    cout<<"Particulas colocadas en la caja: "<<endl;

}


void SetRandomPart2(particula *part, parametros *param)
{
    int nok=0;
    bool good;
    double dis;
    double dmin;
    double diam=1.0;
    double dx,dy,dz;
    //double poly=0.36;
    double Vpart=0.0;
    double vol;
    double u1,u2;
    for(int i=0;i< param->npart;i++)
    {
#ifdef POLYDISPERSE
        u1=double(rand())/RAND_MAX;
        u2=double(rand())/RAND_MAX;
        diam=exp(poly*sqrt(-2.*log(u1))*cos(2.*pi*u2));
#endif
        part[i].diameter=diam;
        part[i].mass=diam*diam*diam;
        Vpart+=1./6.*pi*pow(part[i].diameter,3.0);
    }
    param->vol=Vpart/param->phi;

    param->lbox[0]=pow(param->vol/param->ratio,1./3.);
    param->lbox[1]=param->lbox[0];
    param->lbox[2]=param->ratio*param->lbox[0];
    dmin=param->lbox[0];

    for(int i=0;i< param->npart;i++)
    {
        part[i].pos[0]=param->lbox[0]*(double(rand())/RAND_MAX);
        part[i].pos[1]=param->lbox[1]*(double(rand())/RAND_MAX);
        //part[i].z=bx.lz*(double(rand())/RAND_MAX);
        part[i].pos[2]=(param->lbox[2]-2.*part[i].diameter)*(double(rand())/RAND_MAX)+part[i].diameter;
        //part[i].z=(bx.lz-2.*part[i].sigma)*(double(rand())/RAND_MAX)+part[i].sigma;
        BC(&part[i],*param);
        ///Comprobamos solapamiento con otras particulas
        good=true;
        double daux=1.e7;
        for(int j=0;j<nok;j++)
        {
            ///Recorremos todas las particulas que ya esten metidas en la lista
            if(good)
            {
                dis=distancia(part[i],part[j],&dx,&dy,&dz,*param);
                if(dis<1.1*0.5*(part[i].diameter+part[j].diameter))
                {
                    ///Hay solapamiento entre las particulas i y j. No aceptamos la particula i
                    good=false;
                    i--;
                } 
                else
                {
                        if(dis<daux) daux=dis;
                }
            }
        }
        if(good)
        {
                if( dmin>daux) dmin=daux;
                ///Si pasa por aqui la particula no se solapa con ninguna otra --> se acepta
                nok++;
        }
   }
    cout<<"Particulas colocadas en la caja: (dmin: "<<dmin<<")"<<endl;

}

#endif  /* INICIAL_H */
